clear all; close all; 
cucd=pwd; 
loadSubPath='\BMJ_new\DtrKimbReNormShellReal\pr22 twoWEmp\\end 2007.50 csminwel pr22 TwoWEmp 05 05 14'; 

loadPath=fullfile(cucd,loadSubPath); 

cd(loadPath); 
load workspace; 
cd(cucd); 

Npoints=15; 
perturb.one.name='lambda'; 
perturb.one.grid=linspace(0.67,0.8,Npoints); 

perturb.two.name='STDh'; 
perturb.two.grid=linspace(1.1,1.7,Npoints); 


perturb.one.posFull=cellposition(perturb.one.name,names.param); 
perturb.one.posEst =cellposition(perturb.one.name,names.param(posStru.param.est) ); 

perturb.two.posFull=cellposition(perturb.two.name,names.param); 
perturb.two.posEst =cellposition(perturb.two.name,names.param(posStru.param.est) ); 

%% Display starting values 
perturb.one.estimated=model.param( perturb.one.posFull ); 
perturb.two.estimated=model.param( perturb.two.posFull ); 
dispaj( perturb.one.name,'=',perturb.one.estimated ); 
dispaj( perturb.two.name,'=',perturb.two.estimated ); 

%% Add estimated values to grid 
perturb.one.grid=sort( [perturb.one.grid(:);perturb.one.estimated ] ); 
perturb.two.grid=sort( [perturb.two.grid(:);perturb.two.estimated ] ); 


%% Evaluate Initial Density 
% Extracted the estimate Mode 
modeEstimatedOnly=model.param(posStru.param.est);
% Evaluate the Initial density 
fInitial=...
    feval(filterStru.funcPost,modeEstimatedOnly,parStru,...
    posStru,model,dataStru,prpar,prss,filterStru,flags);
if abs(fInitial-model.logPost) > 1e-4
    warning('Starting Posterior for Hessian and Posterior at Mode do not coincide')
end

modeTemp=modeEstimatedOnly; 

[xmat,ymat]=meshgrid(perturb.one.grid,perturb.two.grid); 
nx=length(xmat); 
ny=length(ymat); 
fmat=nan(nx,ny); 

%%
count=0;
for hh=1:nx*ny 
    
    modeTemp(perturb.one.posEst)=xmat(hh);
    modeTemp(perturb.two.posEst)=ymat(hh);
    dispaj( perturb.one.name,'=',xmat(hh),'  ',perturb.two.name,'=',ymat(hh) ); 
    
    count=count+1;
        fTemp=feval(filterStru.funcPost,modeTemp,parStru,...
    posStru,model,dataStru,prpar,prss,filterStru,flags);
if fTemp > -1e10 
    fmat(hh)=fTemp; 
end 

end

xmat=xmat'; 
ymat=ymat'; 


handle_vec=zeros(2,1); 

handle_vec(1)=figure; 
surf(xmat,ymat,fmat); 
ylabel(perturb.two.name); 
xlabel(perturb.two.name); 

handle_vec(2)=figure; 
CC=contour(xmat,ymat,fmat,15); 
ylabel(perturb.two.name); 
xlabel(perturb.one.name); 
clabel(CC); 